Cargar paquetes

library(readxl)
library(tidyverse)
library(kableExtra)
library(sf)
library(mapview)
library(vegan)
library(devtools)
library(indicspecies)
library(SpadeR)
library(adespatial) #TBI, para comparar matrices de comunidad entre estaciones
library(RColorBrewer)
library(spdep)
download.file('https://raw.githubusercontent.com/biogeografia-master/scripts-de-analisis-BCI/master/biodata/funciones.R', 'funciones1.R', )
source('funciones1.R')
source('funciones2.R')
source('estadistica-zonal.R')

Cargar datos

tabla_odk <- read_csv('tabla ODK.csv')
tabla_odk <- tabla_odk[grep('prueba', tabla_odk$codigo, ignore.case = T, invert = T), ]
tabla_odk$codigo[grep('h14m75', tabla_odk$codigo, ignore.case = T)] <- 'H10M75'
tabla_odk$codigo[grep('h1m14', tabla_odk$codigo, ignore.case = T)] <- 'H12M14'
tabla_odk$codigo <- toupper(tabla_odk$codigo)
tabla_odk$codigo[grep('h1m88', tabla_odk$codigo, ignore.case = T)] <- 'H2M88'
tabla_odk$codigo[grep('2h23h42', tabla_odk$codigo, ignore.case = T)] <- '2H23M42'
tabla_odk$codigo[grep('2h23h43', tabla_odk$codigo, ignore.case = T)] <- '2H23M43'
tabla_odk$codigo[grep('2H25M147', tabla_odk$codigo, ignore.case = T)] <- '2H24M147'

datos_orig <- read_xlsx(
  path = 'Datos de muestreo Marzo y Mayo 2022.xlsx',
  col_types = c('numeric', rep('text', 7), 'date', 'numeric', rep('text', 5)))
which(is.na(match(tabla_odk$codigo, datos_orig$Código))) #Todos los elementos tienen contrapartes entre tabla si "integer(0)"
## integer(0)
ncol(datos_orig) # Número de columnas
## [1] 15
nrow(datos_orig)
## [1] 255
colnames(datos_orig)
##  [1] "Parcela"          "Código"           "Coordenadas"      "Datum"           
##  [5] "Altura"           "Incertidumbre"    "Especie"          "Colector/os"     
##  [9] "Fecha"            "No. de individos" "Observaciones"    "País"            
## [13] "Provincia"        "Municipio"        "Área protegida"
datos_01 <- datos_orig %>%
  mutate(lat = str_split(datos_orig$Coordenadas, pattern = ' ', simplify = T)[,1],
         lon = str_split(datos_orig$Coordenadas, pattern = ' ', simplify = T)[,2]) %>% 
  mutate_at(c('lat', 'lon'), as.numeric)
datos_02 <- datos_01 %>%
  select(-Parcela) %>%
  mutate(Hexágono = as.numeric(gsub('.*H', '', gsub('M.*', '', Código)))) %>% 
  select(Hexágono, everything())
datos_02 %>%
  select(Hexágono, Código, Especie) %>% 
  head() %>% 
  kable()
Hexágono Código Especie
13 H13M1 Pheidole subarmata
13 H13M2 Dorymyrmex antillanus
13 H13M3 Paratrechina longicornis
13 H13M4 Pheidole subarmata
13 H13M5 Solenopsis geminata
13 H13M6 Pheidole jelskii
datos_03 <- datos_02 %>%
  inner_join(
    tabla_odk %>%
      select(codigo, `coordenadaid-Latitude`, `coordenadaid-Longitude`), by = c('Código' = 'codigo')) %>% 
  select(-lat,-lon) %>% 
  rename(lat = `coordenadaid-Latitude`,
         lon = `coordenadaid-Longitude`)
datos_03[grep('^H13M1$', datos_03$Código, ignore.case = T), 'Hexágono'] <- 12 #Misplaced, corrected
datos_03$Código[grep('^H13M1$', datos_03$Código, ignore.case = T)] <- 'H12M1' #Misplaced, corrected
datos_03[grep('^H1M87$', datos_03$Código, ignore.case = T), 'Hexágono'] <- 4 #Misplaced, corrected
datos_03$Código[grep('^H1M87$', datos_03$Código, ignore.case = T)] <- 'H4M87' #Misplaced, corrected
datos_03$Especie <- gsub('Pheidole  jelskii', 'Pheidole jelskii', datos_03$Especie)
datos_03$Especie <- gsub('Ontomachus ruginosa', 'Odontomachus ruginodis', datos_03$Especie)
datos_03$Especie <- gsub('Monomorium floricula', 'Monomorium floricola', datos_03$Especie)
datos_03$Especie <- gsub('Monomorium lanuginosa', 'Tetramorium lanuginosum', datos_03$Especie)
# datos_03 %>% write_csv('datos_03.csv')
# write_excel_csv(datos_03, 'datos_03_compatible_con_excel.csv')

Crear mapa de las colectas

# datos_03 <- read_csv('datos_03.csv')
datos_sf <- datos_03 %>% 
  st_as_sf(coords = c('lon', 'lat'))
st_crs(datos_sf) <- 4326
datos_sf %>% plot()

# st_write(datos_sf, 'datos_sf.gpkg', delete_dsn = T)

# Generar mapa de cuadros sin simbología
## Centroid
# datos_sf %>% mutate(i=1) %>% 
#     group_by(i) %>% 
#     summarize(geometry = st_union(geometry)) %>% 
#     st_centroid %>% mutate(lng )
## Mapa
mapa_cuadros <- mapView(
  datos_sf,
  col.regions = 'grey80',
  alpha.regions = 0.3,
  map.types = 'Esri.WorldImagery',
  legend = F, zoom = 14,
  zcol = 'Código') %>% leafem::addStaticLabels(label = datos_sf$Código) # %>%
  # leaflet::setView(
  #   lng = -79.85136,
  #   lat = 9.15097,
  #   zoom = 15)
mapa_cuadros
# mapa_cuadros %>% mapshot(file = 'mapa_punto.png') #Genera archivo
# Capa de hexagonos
hex <- st_read('cuadricula-final.gpkg')
## Reading layer `cuadricula-final' from data source 
##   `C:\Users\User\Desktop\Universidad\Tesis\Datos de los muestreos\cuadricula-final.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 28 features and 9 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 402697.2 ymin: 2041872 xmax: 403143.5 ymax: 2042260
## Projected CRS: WGS 84 / UTM zone 19N
hex %>% inner_join(#NUMERO DE NIDOS POR HEXAGONO
  datos_03 %>%
    group_by(Hexágono) %>% 
  summarise(n = n()), by = c('ENLACE' = 'Hexágono')) %>% 
  mapview(zcol = 'n')

Tabla de número de nidos por Hexágono

datos_03 %>% group_by(Hexágono) %>%
  summarise(n = n()) %>% arrange(desc(n)) %>% 
  head() %>% 
  kable()
Hexágono n
12 19
16 16
21 16
20 15
13 13
6 12

Crear matriz de comunidad

# datos_03 <- read_csv('datos_03.csv') #Descomentar si se desea comenzar el script desde este punto (no olvidar cargar paquetes)
datos_03
## # A tibble: 255 × 17
##    Hexágono Código Coordenadas            Datum Altura  Incert…¹ Especie Colec…²
##       <dbl> <chr>  <chr>                  <chr> <chr>   <chr>    <chr>   <chr>  
##  1       12 H12M1  18.4675906 -69.9139862 WGS84 53 msnm 5 metros Pheido… D. Guz…
##  2       13 H13M2  18.4672803 -69.9192371 WGS84 54 msnm 4 metros Dorymy… D. Guz…
##  3       13 H13M3  18.4672054 -69.9192181 WGS84 54 msnm 5 metros Paratr… D. Guz…
##  4       13 H13M4  18.4669508 -69.9193136 WGS84 54 msnm 7 metros Pheido… D. Guz…
##  5       13 H13M5  18.4671236 -69.9191236 WGS84 54 msnm 3 metros Soleno… D. Guz…
##  6       13 H13M6  18.4671876 -69.9193009 WGS84 53 msnm 4 metros Pheido… D. Guz…
##  7       12 H12M7  18.4675406 -69.9194121 WGS84 53 msnm 7 metros Pheido… D. Guz…
##  8       12 H12M8  18.4675806 -69.9193971 WGS84 53 msnm 5 metros Pheido… D. Guz…
##  9       12 H12M9  18.4675541 -69.9193724 WGS84 53 msnm 7 metros Pheido… D. Guz…
## 10       12 H12M10 18.4678106 -69.9148187 WGS84 53 msnm 7 metros Soleno… D. Guz…
## # … with 245 more rows, 9 more variables: Fecha <dttm>,
## #   `No. de individos` <dbl>, Observaciones <chr>, País <chr>, Provincia <chr>,
## #   Municipio <chr>, `Área protegida` <chr>, lat <dbl>, lon <dbl>, and
## #   abbreviated variable names ¹​Incertidumbre, ²​`Colector/os`
generar_mc <- function(patron = '^H.*', presencia_ausencia = T){
  datos_03 %>%
    filter(grepl(patron, Código)) %>% 
    select(Hexágono, Especie) %>%
    mutate(n = 1) %>% 
    pivot_wider(names_from = Especie, values_from = n,
                values_fn = ~ sum(.x, na.rm = TRUE),
                values_fill = 0) %>%
    arrange(Hexágono) %>% 
    column_to_rownames('Hexágono') %>% 
    {if(presencia_ausencia) replace(., . > 1, 1) else .}
}
mc_marzo <- generar_mc(patron = '^H.*')
# png('grafico_de_mosaicos_marzo.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(
  add_row(mc_marzo, .before = 1) %>% replace(is.na(.), 0) %>% remove_rownames()
)

# dev.off()
mc_mayo <- generar_mc(patron = '^2H.*')
# png('grafico_de_mosaicos_mayo.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(mc_mayo)

# dev.off()
mc_pooled <- generar_mc(patron = '.*')
# png('grafico_de_mosaicos_pooled.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(mc_pooled)

# dev.off()
mc_marzo_nidos <- generar_mc(patron = '^H.*', presencia_ausencia = F)
# png('grafico_de_mosaicos_nidos_marzo.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(
  add_row(mc_marzo_nidos, .before = 1) %>% replace(is.na(.), 0) %>% remove_rownames()
)

# dev.off()

mc_mayo_nidos <- generar_mc(patron = '^2H.*', presencia_ausencia = F)
# png('grafico_de_mosaicos_nidos_mayo.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(mc_mayo_nidos)

# dev.off()

# Nidos combinados
mc_nidos <- generar_mc(patron = '.*', presencia_ausencia = F)
# png('grafico_de_mosaicos_nidos.png', width = 3840, height = 2160, pointsize = 14, res = 350)
crear_grafico_mosaico_de_mc(mc_nidos)

# dev.off()

# Reescalado
mc_puente <- scales::rescale(as.matrix(mc_nidos), to= c(1, 5))
mc_nidos_reescalado <- round(ifelse(mc_puente==1, 0, mc_puente), 0)

EDA

mc_todas <- list(mc_marzo, mc_mayo, mc_pooled)
names(mc_todas) <- c('Marzo', 'Mayo', 'Combinada')
#' ### Lista de especies
map(mc_todas, 
    ~ sort(colnames(.x)) %>%
      as_tibble() %>%
      mutate(fila = 1:nrow(.)) %>%
      select(2, 1) %>% 
      kable(col.names = c('#', 'Especie'), table.attr='class="myTable"'))
$Marzo
# Especie
1 Cyphomyrmex rimosus
2 Dorymyrmex antillanus
3 Mycetomoellerius jamaicensis
4 Odontomachus ruginodis
5 Paratrechina longicornis
6 Pheidole jelskii
7 Pheidole moerens
8 Pheidole subarmata
9 Solenopsis geminata
10 Tetramorium bicarinatum
11 Tetramorium lanuginosum
12 Wasmannia auropunctata
$Mayo
# Especie
1 Acropyga dubitata
2 Brachymyrmex heeri
3 Cyphomyrmex rimosus
4 Dorymyrmex antillanus
5 Monomorium floricola
6 Mycetomoellerius jamaicensis
7 Odontomachus ruginodis
8 Paratrechina longicornis
9 Pheidole jelskii
10 Pheidole subarmata
11 Pogonomyrmex schmitti
12 Solenopsis geminata
13 Tetramorium lanuginosum
14 Wasmannia auropunctata
$Combinada
# Especie
1 Acropyga dubitata
2 Brachymyrmex heeri
3 Cyphomyrmex rimosus
4 Dorymyrmex antillanus
5 Monomorium floricola
6 Mycetomoellerius jamaicensis
7 Odontomachus ruginodis
8 Paratrechina longicornis
9 Pheidole jelskii
10 Pheidole moerens
11 Pheidole subarmata
12 Pogonomyrmex schmitti
13 Solenopsis geminata
14 Tetramorium bicarinatum
15 Tetramorium lanuginosum
16 Wasmannia auropunctata

Matriz ambiental

cobertura <- st_read('coberturas.gpkg')
## Reading layer `coberturas' from data source 
##   `C:\Users\User\Desktop\Universidad\Tesis\Datos de los muestreos\coberturas.gpkg' 
##   using driver `GPKG'
## Simple feature collection with 82 features and 1 field
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 402699.2 ymin: 2041878 xmax: 403116 ymax: 2042246
## Projected CRS: WGS 84 / UTM zone 19N
hex_inters <- st_as_sf(st_intersection(st_union(cobertura), hex)) %>% mutate(ENLACE = hex$ENLACE)
ma_sf <- zstatsf(zones = hex_inters, values = cobertura, grpx = "ENLACE", grpy = "tipo")
ma_sf <- ma_sf %>% mutate(across(2:5, units::drop_units))
ma_sf %>% plot()

ma <- ma_sf %>% st_drop_geometry() %>% column_to_rownames('ENLACE')
ma
##    construido, mobiliario (bordes edificios, acerado, bancos, postes...)
## 1                                                              0.0000000
## 2                                                             18.4869756
## 3                                                              0.0000000
## 4                                                              3.4888882
## 5                                                              0.0000000
## 6                                                              2.4769313
## 7                                                              0.6085964
## 8                                                             27.2501334
## 9                                                              7.5820449
## 10                                                             0.1769173
## 11                                                             5.1515836
## 12                                                            38.0738217
## 13                                                            51.8993681
## 14                                                             0.0000000
## 15                                                             1.1269871
## 16                                                             3.7630314
## 17                                                            69.1243002
## 18                                                            16.2492695
## 19                                                            25.6929698
## 20                                                            15.8724897
## 21                                                            21.0199767
## 22                                                            43.0273268
## 23                                                            11.3836316
## 24                                                            44.3322007
## 25                                                            23.9599716
## 26                                                             1.7747530
## 27                                                            18.4280883
## 28                                                             0.3534807
##       dosel edificación erguida suelo, herbáceas, no edificado ni cubierto
## 1  89.86997            0.000000                                10.13003022
## 2  80.17277            0.000000                                 1.34025095
## 3  99.88796            0.000000                                 0.11204383
## 4  61.57379           24.838187                                10.09913880
## 5  77.72845           22.271546                                 0.00000000
## 6  95.53464            0.000000                                 1.98843129
## 7  97.26926            0.000000                                 2.12214164
## 8  67.74032            0.000000                                 5.00955070
## 9  54.37792           37.902586                                 0.13744975
## 10 88.68445            0.000000                                11.13863183
## 11 80.03134            1.765534                                13.05153872
## 12 52.72902            0.000000                                 9.19716130
## 13 46.72021            1.068254                                 0.31216408
## 14 96.50561            0.000000                                 3.49438827
## 15 97.93200            0.000000                                 0.94101640
## 16 96.23697            0.000000                                 0.00000000
## 17 26.98363            0.000000                                 3.89206487
## 18 63.99899           19.294636                                 0.45710407
## 19 74.29148            0.000000                                 0.01555076
## 20 84.12751            0.000000                                 0.00000000
## 21 67.26251            0.000000                                11.71751267
## 22 15.30474           28.156893                                13.51103755
## 23 41.44998           46.221870                                 0.94452067
## 24 53.38524            0.000000                                 2.28256075
## 25 76.04003            0.000000                                 0.00000000
## 26 65.53121           24.489490                                 8.20454626
## 27 78.54138            0.000000                                 3.03053153
## 28 93.44066            0.000000                                 6.20586141
rowSums(ma) #Comprobación de suma 100 por filas
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
## 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 
##  21  22  23  24  25  26  27  28 
## 100 100 100 100 100 100 100 100

Matriz combinada

Resúmenes

#' ### Número de sitios, tanto en matriz de comunidad como en ambiental
#' Verifica que coinciden
nrow(mc_todas$Combinada) #En la matriz de comunidad
## [1] 28
nrow(ma) #En la matriz ambiental
## [1] 28
#' ### Riqueza numérica de especies (usando matriz de comunidad) por quadrat
#' Nota: cargar paquete vegan arriba, en el área de paquetes
specnumber(mc_todas$Combinada)
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 
##  2  3  2  4  4  7  9  3  4  2  6  5  6  2  6  3  1  3  4  4  6  3  2  4  5  2 
## 27 28 
##  4  4
sort(specnumber(mc_todas$Combinada)) # Ordenados ascendentemente
## 17  1  3 10 14 23 26  2  8 16 18 22  4  5  9 19 20 24 27 28 12 25 11 13 15 21 
##  1  2  2  2  2  2  2  3  3  3  3  3  4  4  4  4  4  4  4  4  5  5  6  6  6  6 
##  6  7 
##  7  9
summary(specnumber(mc_todas$Combinada)) # Resumen estadístico
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   2.750   4.000   3.929   5.000   9.000
#' ### Riqueza numérica de toda la "comunidad"
specnumber(colSums(mc_todas$Combinada))
## [1] 16

Agrupamiento, asociacion de especies con habitat (Objetivo 3 o 4)

En esta seccion se realiza una clasificacion no supervisada de habitats en funcion de su composicion de especies de hormigas, y se explora cuales atributos (de los 4 generados a partir del mapa de cobertura) caracterizan a dichos habitats.

Filtrado para excluir hexagonos muy pobres o especies muy raras, si aplica

mc_orig <- mc_todas$Combinada
data.frame(`Número de hexágonos` = sort(colSums(mc_orig), decreasing = T), check.names = F) %>% 
  kable(booktabs=T) %>%
  kable_styling(latex_options = c("HOLD_position", "scale_down")) %>%
  gsub(' NA ', '', .) # Número de hexágonos en los que está presente cada especie
Número de hexágonos
Pheidole subarmata 25
Solenopsis geminata 20
Dorymyrmex antillanus 19
Paratrechina longicornis 8
Pheidole jelskii 8
Wasmannia auropunctata 6
Cyphomyrmex rimosus 6
Odontomachus ruginodis 5
Tetramorium lanuginosum 3
Pogonomyrmex schmitti 3
Mycetomoellerius jamaicensis 2
Pheidole moerens 1
Tetramorium bicarinatum 1
Brachymyrmex heeri 1
Monomorium floricola 1
Acropyga dubitata 1
# Usa el vector anterior para determinar un umbral o rango de registros para filtrar tu matriz
# ¿En cuántos hexágonos está cada especie? Filtra tus datos usando tu propio criterio.
# Especies que aparecen en pocos hexágonos se consideran "raras". Por ejemplo, si una especie sólo
# aparece en un hexágono en todo el país, es un "singleton", si en dos, "doubleton", y así.
# Estas especies podrían contribuir a generar "ruido" en análisis posteriores, se recomienda excluirlas.
# Elige un valor mínimo (representado por único número entero) o por un rango de enteros (e.g. de 10 a 20),
# para seleccionar las especies que estén mejor representadas de acuerdo a tu criterio.
# Por ejemplo, si usas el valor m, el script considerará a este valor como "el número mínimo de hexágonos
# en los que está representada una especie, y creará una matriz de comunidad de especies seleccionadas
# que están presentes en m hexágonos o más. Si eliges un rango, por ejemplo [m,n], el script generará
# una matriz de comunidad que representadas un mínimo de m hexágonos y un máximo de n hexágonos.
# (ambos extremos inclusive).
en_cuantos_hex <- 2 #!!!ESTE UMBRAL EXCLUYE A ESPECIES QUE APARECEN EN SOLO 1 HEXAGONO, Y CONSERVA SOLO AQUELLAS QUE APARECEN EN 2 O MAS HEXAGONOS
# Explicación: "en_cuantos_hex <- X", donde X es el número de hexágonos mínimo donde cada especie
# debe estar presente. IMPORTANTE: elige TU PROPIO umbral.
{if(length(en_cuantos_hex)==1) selector <- en_cuantos_hex:max(colSums(mc_orig)) else
  if(length(en_cuantos_hex)==2)
    selector <- min(en_cuantos_hex):max(en_cuantos_hex) else
      stop('Debes indicar uno o dos valores numéricos')}
selector
##  [1]  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
mc_orig_seleccionadas <- mc_orig[, colSums(mc_orig) %in% selector]

# Mínimo número de especies por hexágono
data.frame(`Número de especies por hexágono` = sort(rowSums(mc_orig), decreasing = T), check.names = F) %>% 
  kable(booktabs=T) %>%
  kable_styling(latex_options = c("HOLD_position", "scale_down")) %>%
  gsub(' NA ', '', .) # Número de hexágonos en los que está presente cada especie
Número de especies por hexágono
7 9
6 7
11 6
13 6
15 6
21 6
12 5
25 5
4 4
5 4
9 4
19 4
20 4
24 4
27 4
28 4
2 3
8 3
16 3
18 3
22 3
1 2
3 2
10 2
14 2
23 2
26 2
17 1
min_especies_por_hex <- 2 #!!!ESTA LINEA SACA LOS HEXAGONOS QUE TIENEN SOLO 1 ESPECIE, Y CONSERVA SOLO AQUELLOS QUE TIENEN 2 O MAS ESPECIES PRESENTES (EN EL CASO DE LA MATRIZ DE ESPECIES COMBINADAS, EXCLUYE EL HEXAGONO 17)
# Explicación: "min_especies_por_hex <- Y", donde Y es el número mínimo (inclusive) de especies
# que debe existir en cada hexágono. Por debajo de dicho valor, el hexágono es excluido.
mi_fam <- mc_orig_seleccionadas[rowSums(mc_orig_seleccionadas)>=min_especies_por_hex, ]
nrow(mi_fam)
## [1] 27
# mi_fam <- mc_orig_seleccionadas[!rowSums(mc_orig_seleccionadas)==0, ] #Elimina filas sin registros
# rowSums(mi_fam) #Riqueza por hexágonos con especies seleccionadas. Comentado por extenso
all(rowSums(mi_fam)>0) #Debe ser TRUE: todos los hexágonos tienen al menos 1 registro
## [1] TRUE
ncol(mi_fam) #Riqueza de especies
## [1] 11
# IMPORTANTE, quitar Formicidae. ¡¡Muchos nombres que habría que resolver si se quisiera publicar!!
# Usar nombres cortos o abreviados para las especies
# nombres_largos <- colnames(mi_fam)
# (colnames(mi_fam) <- make.cepnames(word(colnames(mi_fam), 1, 2)))
# (df_equivalencias <- data.frame(
#   nombre_original = nombres_largos,
#   abreviado = colnames(mi_fam)))

Análisis de agrupamiento.

En esta subseccion se realiza propiamente el analisis de agrupamiento segun cuatro metodos distintos (single, completate, upgma, Ward); se elige, al final, el metodo Ward. Se contrastan los metodos con remuestreo multiescalas por bootstrap.

mi_fam_t <- decostand(mi_fam, 'hellinger') #Hellinger
env <- ma[match(rownames(mi_fam), rownames(ma)), ]
all(rownames(mi_fam) == rownames(env)) #Si es TRUE, sigue adelante
## [1] TRUE
library(broom)
library(cluster)
library(gclus)
library(pvclust)
mi_fam_d <- vegdist(mi_fam_t, "euc")
lista_cl <- list(
        cl_single = hclust(mi_fam_d, method = 'single'),
        cl_complete = hclust(mi_fam_d, method = 'complete'),
        cl_upgma = hclust(mi_fam_d, method = 'average'),
        cl_ward = hclust(mi_fam_d, method = 'ward.D2')
)
par(mfrow = c(2,2))
invisible(map(names(lista_cl), function(x) plot(lista_cl[[x]], main = x, hang = -1, cex = 1)))

par(mfrow = c(1,1))
cutree(lista_cl$cl_ward, 2)
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 18 19 20 21 22 23 24 25 26 27 
##  1  1  2  1  1  2  2  2  2  1  2  1  1  1  2  1  1  1  1  1  1  1  1  2  2  1 
## 28 
##  2

Remuestreo multiescalar por bootstrap

if(interactive()) dev.new()
cl_pvclust_ward <-
        pvclust(t(mi_fam_t),
                method.hclust = "ward.D2",
                method.dist = "euc",
                iseed = 999, # Resultado reproducible
                parallel = TRUE)
## Creating a temporary cluster...done:
## socket cluster with 3 nodes on host 'localhost'
## Multiscale bootstrap... Done.
# Añadir los valores de p
plot(cl_pvclust_ward, hang = -1)
# Añadir rectángulos a los grupos significativos
lines(cl_pvclust_ward)
pvrect(cl_pvclust_ward, alpha = 0.91, border = 4)

Guardando agrupamientos a archivos

w_dend_reord <- reorder.hclust(lista_cl$cl_ward, mi_fam_d)
grupos_ward_k2 <- as.factor(cutree(lista_cl$cl_ward, k = 2))
table(grupos_ward_k2) #¿Cuántos hexágonos hay en cada grupo?
## grupos_ward_k2
##  1  2 
## 17 10
plot(w_dend_reord, hang = -1)
rect.hclust(tree = w_dend_reord, k = 2)

grupos_ward_k3 <- as.factor(cutree(lista_cl$cl_ward, k = 3))
table(grupos_ward_k3) #¿Cuántos hexágonos hay en cada grupo?
## grupos_ward_k3
##  1  2  3 
##  6 10 11
plot(w_dend_reord, hang = -1)
rect.hclust(tree = w_dend_reord, k = 3)

# Guardaré estos vectores en archivos para reutilizarlos en *scripts* posteriores: 
# saveRDS(grupos_ward_k2, 'grupos_ward_k2.RDS')
# saveRDS(grupos_ward_k3, 'grupos_ward_k3.RDS') #No hubo patrones significativos de asociacion

Caractersticas ambientales de grupos, asociacion con habitat

ma_para_combinadas <- ma[match(rownames(mi_fam), rownames(ma)), ]
(m_amb_ward_k2 <- ma_para_combinadas %>%
   rownames_to_column('hex_id') %>% 
   mutate(grupos_ward_k2) %>%
   pivot_longer(-c(grupos_ward_k2, hex_id), names_to = "variable", values_to = "valor"))
## # A tibble: 108 × 4
##    hex_id grupos_ward_k2 variable                                          valor
##    <chr>  <fct>          <chr>                                             <dbl>
##  1 1      1              construido, mobiliario (bordes edificios, acerad…  0   
##  2 1      1              dosel                                             89.9 
##  3 1      1              edificación erguida                                0   
##  4 1      1              suelo, herbáceas, no edificado ni cubierto        10.1 
##  5 2      1              construido, mobiliario (bordes edificios, acerad… 18.5 
##  6 2      1              dosel                                             80.2 
##  7 2      1              edificación erguida                                0   
##  8 2      1              suelo, herbáceas, no edificado ni cubierto         1.34
##  9 3      2              construido, mobiliario (bordes edificios, acerad…  0   
## 10 3      2              dosel                                             99.9 
## # … with 98 more rows
m_amb_ward_k2_tw <- m_amb_ward_k2 %>%
  group_by(variable) %>%
  summarise(
    p_valor_t = tryCatch(t.test(valor ~ grupos_ward_k2)$p.value, error = function(e) NA),
    p_valor_w = tryCatch(wilcox.test(valor ~ grupos_ward_k2)$p.value, error = function(e) NA)
    ) %>%
  arrange(p_valor_t)
m_amb_ward_k2_tw %>%
  kable(booktabs=T) %>%
  kable_styling(latex_options = c("HOLD_position", "scale_down")) %>%
  gsub(' NA ', '', .)
variable p_valor_t p_valor_w
construido, mobiliario (bordes edificios, acerado, bancos, postes…) 0.0405413 0.1910607
dosel 0.0696707 0.0743239
suelo, herbáceas, no edificado ni cubierto 0.6704152 0.8405846
edificación erguida 0.7283481 0.8108505
m_amb_ward_k2 %>% 
  group_by(variable) %>% 
        ggplot() + aes(x = grupos_ward_k2, y = valor, fill = grupos_ward_k2) + 
        geom_boxplot(lwd = 0.2) + 
        scale_fill_brewer(palette = 'Set1') +
        theme_bw(base_size=12) +
        theme(legend.position="none") +
        facet_wrap(~ variable, scales = 'free_y', ncol = 2,
                   labeller = label_wrap_gen(width = 35,multi_line = T))

# CON WARD K3 NO APARECEN EFECTOS SIGNIFICATIVOS. MEJOR USAR SOLO AGRUPAMIENTO WARD K2

En resumen, las variables de porcentaje de “construido…” y “dosel…” son las unicas que presentaron efectos significativos entre grupos, por lo que consideramos que estas son las que, con mayor probabilidad, se asocian con la composicion de especies diferenciada entre grupos. El grupo 1 se caracteriza por tener areas construidas predominantemente, y el 2 por tener predominio de dosel.

Realizar analisis de asociacion propiamente (correlacion punto biserial, funcion r.g)

phi_ward_k2 <- multipatt(
  mi_fam,
  grupos_ward_k2,
  func = "r.g",
  max.order = 1,
  control = how(nperm = 999))
summary(phi_ward_k2)
## 
##  Multilevel pattern analysis
##  ---------------------------
## 
##  Association function: r.g
##  Significance level (alpha): 0.05
## 
##  Total number of species: 11
##  Selected number of species: 3 
##  Number of species associated to 1 group: 3 
## 
##  List of species associated to each combination: 
## 
##  Group 1  #sps.  1 
##                        stat p.value    
## Dorymyrmex antillanus 0.749   0.001 ***
## 
##  Group 2  #sps.  2 
##                         stat p.value   
## Wasmannia auropunctata 0.655   0.002 **
## Cyphomyrmex rimosus    0.492   0.016 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ma_sf_grupos_ward_k2 <- ma_sf %>%
  inner_join(m_amb_ward_k2 %>%
               select(ENLACE=hex_id, grupos_ward_k2) %>%
               distinct() %>% 
               mutate(ENLACE=as.integer(ENLACE)))
st_geometry(ma_sf_grupos_ward_k2) <- 'geometry'
mapa_ma_sf_grupos_ward_k2 <- ma_sf_grupos_ward_k2 %>% 
  ggplot + aes(fill = grupos_ward_k2) +
  geom_sf()
if(interactive()) dev.new()
mapa_ma_sf_grupos_ward_k2

mapa_ma_sf <- ma_sf_grupos_ward_k2 %>% 
  select(-grupos_ward_k2) %>% 
  pivot_longer(-c(geometry, ENLACE)) %>% 
  group_by(name) %>% 
  mutate(value=scale(value)) %>% 
  ggplot + aes(fill = value) +
  scale_fill_viridis_b() +
  geom_sf() +
  facet_wrap(~name, ncol = 2)
if(interactive()) dev.new()
mapa_ma_sf

Análisis de ordenación

Analisis de ordenacion. No aporta patrones de interes, puesto que la inercia restringida es muy pequena (es habitual cuando se trabaja con datos de presencia/ausencia). Analisis descartado en este caso.

# mi_fam_t_rda <- rda(mi_fam_t ~ .,
#                     ma[-17,] %>% select(all_of(matches('construido|dosel'))),
#                     scale = T)
# summary(mi_fam_t_rda)
# RsquareAdj(mi_fam_t_rda)$adj.r.squared
# vif.cca(mi_fam_t_rda)
# escalado <- 1
# plot(mi_fam_t_rda,
#      scaling = escalado,
#      display = c("sp", "lc", "cn"),
#      main = paste("Triplot de RDA especies ~ var. GSL + ESA, escalamiento", escalado)
# )
# mi_fam_t_rda_sc1 <- scores(mi_fam_t_rda,
#          choices = 1:2,
#          scaling = escalado,
#          display = "sp"
#   )
# # text(mi_fam_t_rda, "species", col="red", cex=0.8, scaling=escalado)
# arrows(0, 0,
#        mi_fam_t_rda_sc1[, 1] * 0.9,
#        mi_fam_t_rda_sc1[, 2] * 0.9,
#        length = 0,
#        lty = 1,
#        col = "red"
# )

Diversidad, variabilidad estacional

Estimacion de riqueza (Chao)

Matriz combinada, ambas estaciones

Pool de especies, estimadores tradicionales de la riqueza (Chao, Jack1, Jack2, etc.). Esta matriz tiene muchas especies raras, por lo que los estimadores penalizan la riqueza capturada. Se considerara usar una matriz ponderada por el numero de nidos.

# mc_orig <- mc_todas$Combinada
specpool(mc_orig)
##     Species     chao  chao.se    jack1 jack1.se    jack2    boot  boot.se  n
## All      16 28.05357 16.54004 20.82143 2.564723 24.57011 18.0219 1.375846 28
specpool(mc_orig)[[1]]/specpool(mc_orig)[-c(3,5,8)]*100
##     Species     chao    jack1    jack2     boot        n
## All     100 57.03374 76.84391 65.11978 88.78088 57.14286
# Nidos
estimateR(colSums(mc_nidos))
##     S.obs   S.chao1  se.chao1     S.ACE    se.ACE 
## 16.000000 26.000000 10.258674 20.465909  2.201645
estimateR(colSums(mc_nidos))[[1]]/estimateR(colSums(mc_nidos))[-c(3,5,8)]*100
##     S.obs   S.chao1     S.ACE 
## 100.00000  61.53846  78.17879

Estimadores

Matriz original combinada: el error “Error in if (var_iChao2 > 0) { : missing value where TRUE/FALSE needed” impide ejecutar la estimacion para la matriz combinada.

# Incidencia
ChaoSpecies(data.frame(V1 = c(nrow(mc_orig), as.integer(colSums(mc_orig)))),
            datatype = 'incidence_freq', k=1, conf=0.95)
# Nmero de nidos
ChaoSpecies(data.frame(V1 = c(nrow(mc_nidos), as.integer(colSums(mc_nidos)))),
            datatype = 'abundance', k=1, conf=0.95)
## Warning: In this case, it can't estimate the variance of Homogeneous estimation
## 
## (1) BASIC DATA INFORMATION:
## 
##                                          Variable Value
##     Sample size                                 n   283
##     Number of observed species                  D    17
##     Coverage estimate for entire dataset        C 0.982
##     CV for entire dataset                      CV 1.538
##     Cut-off point                               k     1
## 
##                                                       Variable Value
##     Number of observed individuals for rare group       n_rare     5
##     Number of observed species for rare group           D_rare     5
##     Estimate of the sample coverage for rare group      C_rare     0
##     Estimate of CV for rare group in ACE               CV_rare     0
##     Estimate of CV1 for rare group in ACE-1           CV1_rare     0
##     Number of observed individuals for abundant group   n_abun   278
##     Number of observed species for abundant group       D_abun    12
## 
## NULL
## 
## 
## (2) SPECIES RICHNESS ESTIMATORS TABLE:
## 
##                               Estimate   s.e. 95%Lower 95%Upper
##     Homogeneous Model              Inf     NA       NA       NA
##     Homogeneous (MLE)           17.000  1.123   17.000   21.705
##     Chao1 (Chao, 1984)          26.965 10.235   18.888   69.583
##     Chao1-bc                    26.965 10.235   18.888   69.583
##     iChao1 (Chiu et al. 2014)   29.465 10.607   19.934   69.955
##     ACE (Chao & Lee, 1992)      26.965 10.235   18.888   69.583
##     ACE-1 (Chao & Lee, 1992)    26.965 10.235   18.888   69.583
##     1st order jackknife         21.982  3.154   18.597   32.544
##     2nd order jackknife         26.947  5.453   20.642   44.166
## 
## 
## (3) DESCRIPTION OF ESTIMATORS/MODELS:
## 
## Homogeneous Model: This model assumes that all species have the same incidence or detection probabilities. See Eq. (3.2) of Lee and Chao (1994) or Eq. (12a) in Chao and Chiu (2016b).
## 
## Chao2 (Chao, 1987): This approach uses the frequencies of uniques and duplicates to estimate the number of undetected species; see Chao (1987) or Eq. (11a) in Chao and Chiu (2016b).
##      
## Chao2-bc: A bias-corrected form for the Chao2 estimator; see Chao (2005).
##   
## iChao2: An improved Chao2 estimator; see Chiu et al. (2014).
## 
## ICE (Incidence-based Coverage Estimator): A non-parametric estimator originally proposed by Lee and Chao (1994) in the context of capture-recapture data analysis. The observed species are separated as frequent and infrequent species groups; only data in the infrequent group are used to estimate the number of undetected species. The estimated CV for species in the infrequent group characterizes the degree of heterogeneity among species incidence probabilities. See Eq. (12b) of Chao and Chiu (2016b), which is an improved version of Eq. (3.18) in Lee and Chao (1994). This model is also called Model(h) in capture-recapture literature where h denotes "heterogeneity".
## 
## ICE-1: A modified ICE for highly-heterogeneous cases.
## 
## 1st order jackknife: It uses the frequency of uniques to estimate the number of undetected species; see Burnham and Overton (1978).
## 
## 2nd order jackknife: It uses the frequencies of uniques and duplicates to estimate the number of undetected species; see Burnham and Overton (1978).
## 
## 95% Confidence interval: A log-transformation is used for all estimators so that the lower bound of the resulting interval is at least the number of observed species. See Chao (1987).

Usando matriz combinada, incluyendo singletons

Para fines didacticos, representamos la curva de acumulacion de especies de la matriz combinada, donde se observa una pronunciada pendiente y una escasa capacidad de alcanzarse la asintota. Decimos “para fines didacticos”, porque complementa el analisis posterior por grupos segun Ward.

mc_orig_lista <- sapply(rownames(mc_orig), function(x) as.numeric(mc_orig[x,]), simplify = F)
nasin_raref_pooled <- iNEXT::iNEXT(
  x = colSums(mc_orig), #if(is.list(mc_orig)) mc_orig_lista else mc_orig,
  q=0,
  knots = 400,
  datatype = 'incidence_freq')
acumulacion_especies_pooled <- iNEXT::ggiNEXT(nasin_raref_pooled, type=1) +
  theme_bw() +
  theme(
    text = element_text(size = 20),
    panel.background = element_rect(fill = 'white', colour = 'black'),
    panel.grid.major = element_line(colour = "grey", linetype = "dashed", size = 0.25)
  ) +
  ylab('Riqueza de especies') +
  xlab('Número de sitios') +
  scale_y_continuous(breaks = seq(0,80, length.out = 9)) +
  scale_color_manual(values = brewer.pal(8, 'Set2')) +
  scale_fill_manual(values = brewer.pal(8, 'Set2'))
if(interactive()) dev.new()
acumulacion_especies_pooled

Numero de nidos

mc_nidos_lista <- sapply(rownames(mc_nidos), function(x) as.numeric(mc_nidos[x,]), simplify = F)
nasin_raref_nidos <- iNEXT::iNEXT(
  x = colSums(mc_nidos), #if(is.list(mc_nidos)) mc_nidos_lista else mc_nidos,
  q=0,
  knots = 400,
  datatype = 'abundance')
acumulacion_especies_nidos <- iNEXT::ggiNEXT(nasin_raref_nidos, type=1) +
  theme_bw() +
  theme(
    text = element_text(size = 20),
    panel.background = element_rect(fill = 'white', colour = 'black'),
    panel.grid.major = element_line(colour = "grey", linetype = "dashed", size = 0.25)
  ) +
  ylab('Riqueza de especies') +
  xlab('Número de sitios') +
  scale_y_continuous(breaks = seq(0,80, length.out = 9)) +
  scale_color_manual(values = brewer.pal(8, 'Set2')) +
  scale_fill_manual(values = brewer.pal(8, 'Set2'))
if(interactive()) dev.new()
acumulacion_especies_nidos

Numero de nidos reescalado

mc_nidos_reescalado_lista <- sapply(rownames(mc_nidos_reescalado), function(x) as.numeric(mc_nidos_reescalado[x,]), simplify = F)
nasin_raref_nidos_reescalado <- iNEXT::iNEXT(
  x = colSums(mc_nidos_reescalado), #if(is.list(mc_nidos)) mc_nidos_lista else mc_nidos,
  q=0,
  knots = 400,
  datatype = 'abundance')
acumulacion_especies_nidos_reescalado <- iNEXT::ggiNEXT(nasin_raref_nidos_reescalado, type=1) +
  theme_bw() +
  theme(
    text = element_text(size = 20),
    panel.background = element_rect(fill = 'white', colour = 'black'),
    panel.grid.major = element_line(colour = "grey", linetype = "dashed", size = 0.25)
  ) +
  ylab('Riqueza de especies') +
  xlab('Número de sitios') +
  scale_y_continuous(breaks = seq(0,80, length.out = 9)) +
  scale_color_manual(values = brewer.pal(8, 'Set2')) +
  scale_fill_manual(values = brewer.pal(8, 'Set2'))
if(interactive()) dev.new()
acumulacion_especies_nidos_reescalado

Usando matriz sin especies singletons, usada en analisis de agrupamiento

Se refiere a la matriz filtrada donde se eliminaron especies poco representadas

specpool(mi_fam)
##     Species chao chao.se jack1 jack1.se    jack2     boot   boot.se  n
## All      11   11       0    11        0 10.10969 11.21473 0.4439909 27
ChaoSpecies(data.frame(V1 = c(nrow(mi_fam), as.numeric(colSums(mi_fam)))),
            datatype = 'incidence_freq', k=2, conf=0.95)
## Warning: In this case, it can't estimate the variance of Homogeneous estimation 
## 
## Warning: In this case, it can't estimate the variance of Model(h) estimation 
## 
## Warning: In this case, it can't estimate the variance of Model(h)-1 estimation 
## 
## Warning: In this case, it can't estimate the variance of 1st-order-jackknife estimation 
## 
## Warning: In this case, it can't estimate the variance of 2nd-order-jackknife estimation 
## 
## Warning: In this case, it can't estimate the variance of Model(h) estimation 
## 
## Warning: In this case, it can't estimate the variance of Model(h)-1 estimation
## 
## (1) BASIC DATA INFORMATION:
## 
##                                          Variable Value
##     Number of observed species                  D    11
##     Number of sampling units                    T    27
##     Total number of incidences                  U   104
##     Coverage estimate for entire dataset        C     1
##     CV for entire dataset                      CV 0.759
## 
##                                                       Variable Value
##     Cut-off point                                            k     2
##     Total number of incidences in infrequent group    U_infreq     2
##     Number of observed species for infrequent group   D_infreq     1
##     Estimated sample coverage for infrequent group    C_infreq     1
##     Estimated CV for infrequent group in ICE         CV_infreq 0.196
##     Estimated CV1 for infrequent group in ICE-1     CV1_infreq 0.196
##     Number of observed species for frequent group       D_freq    10
## 
##                            Q1 Q2
##     Incidence freq. counts  0  1
## 
## 
## (2) SPECIES RICHNESS ESTIMATORS TABLE:
## 
##                               Estimate  s.e. 95%Lower 95%Upper
##     Homogeneous Model            11.00 0.457       11   12.169
##     Chao2 (Chao, 1987)           11.00 0.457       11   12.169
##     Chao2-bc                     11.00 0.457       11   12.169
##     iChao2 (Chiu et al. 2014)    11.00 0.457       11   12.169
##     ICE (Lee & Chao, 1994)       11.00 0.457       11   12.169
##     ICE-1 (Lee & Chao, 1994)     11.00 0.457       11   12.169
##     1st order jackknife          11.00 0.457       11   12.169
##     2nd order jackknife          10.11    NA       NA       NA
## 
## 
## (3) DESCRIPTION OF ESTIMATORS/MODELS:
## 
## Homogeneous Model: This model assumes that all species have the same incidence or detection probabilities. See Eq. (3.2) of Lee and Chao (1994) or Eq. (12a) in Chao and Chiu (2016b).
## 
## Chao2 (Chao, 1987): This approach uses the frequencies of uniques and duplicates to estimate the number of undetected species; see Chao (1987) or Eq. (11a) in Chao and Chiu (2016b).
##      
## Chao2-bc: A bias-corrected form for the Chao2 estimator; see Chao (2005).
##   
## iChao2: An improved Chao2 estimator; see Chiu et al. (2014).
## 
## ICE (Incidence-based Coverage Estimator): A non-parametric estimator originally proposed by Lee and Chao (1994) in the context of capture-recapture data analysis. The observed species are separated as frequent and infrequent species groups; only data in the infrequent group are used to estimate the number of undetected species. The estimated CV for species in the infrequent group characterizes the degree of heterogeneity among species incidence probabilities. See Eq. (12b) of Chao and Chiu (2016b), which is an improved version of Eq. (3.18) in Lee and Chao (1994). This model is also called Model(h) in capture-recapture literature where h denotes "heterogeneity".
## 
## ICE-1: A modified ICE for highly-heterogeneous cases.
## 
## 1st order jackknife: It uses the frequency of uniques to estimate the number of undetected species; see Burnham and Overton (1978).
## 
## 2nd order jackknife: It uses the frequencies of uniques and duplicates to estimate the number of undetected species; see Burnham and Overton (1978).
## 
## 95% Confidence interval: A log-transformation is used for all estimators so that the lower bound of the resulting interval is at least the number of observed species. See Chao (1987).

Curva de acumulacion de especies. Primero generamos la matriz combinada (pooled) por grupos del agrupamiento Ward K2.

mi_fam_ward <- mi_fam %>%
  mutate(g = grupos_ward_k2) %>%
  group_by(g) %>%
  summarise_all(sum) %>%
  select(-g) %>% 
  mutate(N = nrow(mi_fam)) %>% 
  relocate(N, .before = 1) %>% 
  data.frame
mi_fam_ward
##    N Pheidole.subarmata Dorymyrmex.antillanus Paratrechina.longicornis
## 1 27                 16                    16                        3
## 2 27                  9                     2                        5
##   Solenopsis.geminata Pheidole.jelskii Mycetomoellerius.jamaicensis
## 1                  11                6                            1
## 2                   9                2                            1
##   Tetramorium.lanuginosum Odontomachus.ruginodis Wasmannia.auropunctata
## 1                       1                      2                      0
## 2                       2                      3                      6
##   Cyphomyrmex.rimosus Pogonomyrmex.schmitti
## 1                   1                     1
## 2                   5                     2

Posteriormente, generamos la curva de acumlacion.

nasin_raref <- iNEXT::iNEXT(
  x = t(mi_fam_ward),
  q=0,
  knots = 400,
  datatype = 'incidence_freq')
acumulacion_especies <- iNEXT::ggiNEXT(nasin_raref, type=1) +
  theme_bw() +
  theme(
    text = element_text(size = 20),
    panel.background = element_rect(fill = 'white', colour = 'black'),
    panel.grid.major = element_line(colour = "grey", linetype = "dashed", size = 0.25)
  ) +
  ylab('Riqueza de especies') +
  xlab('Número de sitios') +
  scale_y_continuous(breaks = seq(0,80, length.out = 9)) +
  scale_color_manual(values = brewer.pal(8, 'Set2')) +
  scale_fill_manual(values = brewer.pal(8, 'Set2'))
if(interactive()) dev.new()
acumulacion_especies

Por estaciones

Temporada seca

#' ### Número de sitios, tanto en matriz de comunidad como en ambiental
#' Verifica que coinciden
nrow(mc_todas$Marzo) #En la matriz de comunidad
## [1] 27
nrow(ma) #En la matriz ambiental
## [1] 28
#' ### Riqueza numérica de especies (usando matriz de comunidad) por quadrat
#' Nota: cargar paquete vegan arriba, en el área de paquetes
specnumber(mc_todas$Marzo)
##  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 
##  2  1  3  3  2  2  2  2  1  2  5  5  1  2  2  1  3  3  3  4  2  1  2  4  2  2 
## 28 
##  1
sort(specnumber(mc_todas$Marzo)) # Ordenados ascendentemente
##  3 10 14 17 23 28  2  6  7  8  9 11 15 16 22 24 26 27  4  5 18 19 20 21 25 12 
##  1  1  1  1  1  1  2  2  2  2  2  2  2  2  2  2  2  2  3  3  3  3  3  4  4  5 
## 13 
##  5
summary(specnumber(mc_todas$Marzo)) # Resumen estadístico
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.000   2.000   2.000   2.333   3.000   5.000
#' ### Riqueza numérica de toda la "comunidad"
specnumber(colSums(mc_todas$Marzo))
## [1] 12

Comparar matrices entre estaciones con TBI

# TBI
mc_marzo_para_tbi <- cbind(
  rbind(rownames(mc_pooled)[!rownames(mc_pooled) %in% rownames(mc_marzo)], mc_marzo) %>% sapply(., as.integer),
  mc_mayo[,!colnames(mc_mayo) %in% colnames(mc_marzo)] %>% mutate_all(.funs = ~ifelse(.>0, 0, 0)))
colnames(mc_marzo)[!colnames(mc_marzo) %in% colnames(mc_mayo)] #Presentes en marzo, ausentes en mayo
## [1] "Pheidole moerens"        "Tetramorium bicarinatum"
colnames(mc_mayo)[!colnames(mc_mayo) %in% colnames(mc_marzo)] #Presentes en mayo, ausentes en marzo
## [1] "Brachymyrmex heeri"    "Pogonomyrmex schmitti" "Monomorium floricola" 
## [4] "Acropyga dubitata"
mc_mayo_para_tbi <- cbind(
  mc_mayo,
  mc_marzo_para_tbi[,!colnames(mc_marzo_para_tbi) %in% colnames(mc_mayo)] %>% mutate_all(.funs = ~ifelse(.>0, 0, 0)))
mc_mayo_para_tbi <- mc_mayo_para_tbi[,colnames(mc_marzo_para_tbi)]
all.equal(colnames(mc_marzo_para_tbi), colnames(mc_mayo_para_tbi))
## [1] TRUE
rownames(mc_marzo_para_tbi) == rownames(mc_mayo_para_tbi)
##  [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
tbi_porc_diferencias <- TBI(mc_marzo_para_tbi, mc_mayo_para_tbi, method = '%difference', nperm = 999, test.t.perm = T)
tbi_porc_diferencias
## $TBI
##  [1] 0.7142857 0.5000000 1.0000000 0.3333333 0.3333333 0.5555556 1.0000000
##  [8] 0.5000000 0.6000000 0.3333333 0.5000000 0.4285714 0.3333333 1.0000000
## [15] 0.5000000 0.5000000 0.0000000 0.2000000 0.6000000 0.3333333 0.3333333
## [22] 0.2000000 1.0000000 0.6000000 0.4285714 0.0000000 0.3333333 0.6000000
## 
## $p.TBI
##  [1] 0.244 0.592 0.153 0.820 0.803 0.428 0.175 0.576 0.404 0.828 0.574 0.671
## [13] 0.820 0.171 0.568 0.538 0.997 0.912 0.412 0.828 0.823 0.931 0.173 0.382
## [25] 0.645 1.000 0.805 0.397
## 
## $p.adj
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 
## $BCD.mat
##         B/(2A+B+C) C/(2A+B+C) D=(B+C)/(2A+B+C) Change
## Site.1   0.7142857  0.0000000        0.7142857    -  
## Site.2   0.2500000  0.2500000        0.5000000    0  
## Site.3   0.5000000  0.5000000        1.0000000    0  
## Site.4   0.1666667  0.1666667        0.3333333    0  
## Site.5   0.1666667  0.1666667        0.3333333    0  
## Site.6   0.0000000  0.5555556        0.5555556    +  
## Site.7   0.2222222  0.7777778        1.0000000    +  
## Site.8   0.2500000  0.2500000        0.5000000    0  
## Site.9   0.2000000  0.4000000        0.6000000    +  
## Site.10  0.0000000  0.3333333        0.3333333    +  
## Site.11  0.0000000  0.5000000        0.5000000    +  
## Site.12  0.4285714  0.0000000        0.4285714    -  
## Site.13  0.2222222  0.1111111        0.3333333    -  
## Site.14  0.5000000  0.5000000        1.0000000    0  
## Site.15  0.0000000  0.5000000        0.5000000    +  
## Site.16  0.2500000  0.2500000        0.5000000    0  
## Site.17  0.0000000  0.0000000        0.0000000    0  
## Site.18  0.2000000  0.0000000        0.2000000    -  
## Site.19  0.4000000  0.2000000        0.6000000    -  
## Site.20  0.1666667  0.1666667        0.3333333    0  
## Site.21  0.1111111  0.2222222        0.3333333    +  
## Site.22  0.0000000  0.2000000        0.2000000    +  
## Site.23  0.5000000  0.5000000        1.0000000    0  
## Site.24  0.2000000  0.4000000        0.6000000    +  
## Site.25  0.2857143  0.1428571        0.4285714    -  
## Site.26  0.0000000  0.0000000        0.0000000    0  
## Site.27  0.0000000  0.3333333        0.3333333    +  
## Site.28  0.0000000  0.6000000        0.6000000    +  
## 
## $BCD.summary
##  mean(B/den) mean(C/den)   mean(D)   B/(B+C)   C/(B+C) Change
##    0.2047902   0.2866497 0.4914399 0.4167147 0.5832853    +  
## 
## $t.test_B.C
##                 mean(C-B)     Stat   p.param p.perm   p<=0.05
## Paired t.test  0.08185941 1.437708 0.1620069  0.174          
## 
## $BC
## [1] NA
## 
## attr(,"class")
## [1] "TBI"
tbi_jaccard <- TBI(mc_marzo_para_tbi, mc_mayo_para_tbi, method = 'jaccard', nperm = 999, test.t.perm = T)
tbi_jaccard
## $TBI
##  [1] 0.8333333 0.6666667 1.0000000 0.5000000 0.5000000 0.7142857 1.0000000
##  [8] 0.6666667 0.7500000 0.5000000 0.6666667 0.6000000 0.5000000 1.0000000
## [15] 0.6666667 0.6666667 0.0000000 0.3333333 0.7500000 0.5000000 0.5000000
## [22] 0.3333333 1.0000000 0.7500000 0.6000000 0.0000000 0.5000000 0.7500000
## 
## $p.TBI
##  [1] 0.210 0.582 0.171 0.818 0.806 0.395 0.166 0.553 0.397 0.820 0.580 0.656
## [13] 0.795 0.170 0.556 0.554 0.999 0.919 0.435 0.817 0.829 0.923 0.163 0.398
## [25] 0.644 0.999 0.825 0.400
## 
## $p.adj
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 
## $BCD.mat
##         B/(A+B+C) C/(A+B+C) D=(B+C)/(A+B+C) Change
## Site.1  0.8333333 0.0000000       0.8333333    -  
## Site.2  0.3333333 0.3333333       0.6666667    0  
## Site.3  0.5000000 0.5000000       1.0000000    0  
## Site.4  0.2500000 0.2500000       0.5000000    0  
## Site.5  0.2500000 0.2500000       0.5000000    0  
## Site.6  0.0000000 0.7142857       0.7142857    +  
## Site.7  0.2222222 0.7777778       1.0000000    +  
## Site.8  0.3333333 0.3333333       0.6666667    0  
## Site.9  0.2500000 0.5000000       0.7500000    +  
## Site.10 0.0000000 0.5000000       0.5000000    +  
## Site.11 0.0000000 0.6666667       0.6666667    +  
## Site.12 0.6000000 0.0000000       0.6000000    -  
## Site.13 0.3333333 0.1666667       0.5000000    -  
## Site.14 0.5000000 0.5000000       1.0000000    0  
## Site.15 0.0000000 0.6666667       0.6666667    +  
## Site.16 0.3333333 0.3333333       0.6666667    0  
## Site.17 0.0000000 0.0000000       0.0000000    0  
## Site.18 0.3333333 0.0000000       0.3333333    -  
## Site.19 0.5000000 0.2500000       0.7500000    -  
## Site.20 0.2500000 0.2500000       0.5000000    0  
## Site.21 0.1666667 0.3333333       0.5000000    +  
## Site.22 0.0000000 0.3333333       0.3333333    +  
## Site.23 0.5000000 0.5000000       1.0000000    0  
## Site.24 0.2500000 0.5000000       0.7500000    +  
## Site.25 0.4000000 0.2000000       0.6000000    -  
## Site.26 0.0000000 0.0000000       0.0000000    0  
## Site.27 0.0000000 0.5000000       0.5000000    +  
## Site.28 0.0000000 0.7500000       0.7500000    +  
## 
## $BCD.summary
##  mean(B/den) mean(C/den)   mean(D)   B/(B+C)   C/(B+C) Change
##    0.2549603   0.3610261 0.6159864 0.4139058 0.5860942    +  
## 
## $t.test_B.C
##                 mean(C-B)     Stat   p.param p.perm   p<=0.05
## Paired t.test   0.1060658 1.455444 0.1570758  0.152          
## 
## $BC
## [1] NA
## 
## attr(,"class")
## [1] "TBI"
tbi_sorensen <- TBI(mc_marzo_para_tbi, mc_mayo_para_tbi, method = 'sorensen', nperm = 999, test.t.perm = T)
tbi_sorensen
## $TBI
##  [1] 0.7142857 0.5000000 1.0000000 0.3333333 0.3333333 0.5555556 1.0000000
##  [8] 0.5000000 0.6000000 0.3333333 0.5000000 0.4285714 0.3333333 1.0000000
## [15] 0.5000000 0.5000000 0.0000000 0.2000000 0.6000000 0.3333333 0.3333333
## [22] 0.2000000 1.0000000 0.6000000 0.4285714 0.0000000 0.3333333 0.6000000
## 
## $p.TBI
##  [1] 0.219 0.580 0.170 0.805 0.814 0.420 0.157 0.564 0.415 0.819 0.554 0.643
## [13] 0.814 0.167 0.560 0.565 1.000 0.940 0.435 0.827 0.830 0.923 0.148 0.410
## [25] 0.622 1.000 0.789 0.400
## 
## $p.adj
##  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## 
## $BCD.mat
##         B/(2A+B+C) C/(2A+B+C) D=(B+C)/(2A+B+C) Change
## Site.1   0.7142857  0.0000000        0.7142857    -  
## Site.2   0.2500000  0.2500000        0.5000000    0  
## Site.3   0.5000000  0.5000000        1.0000000    0  
## Site.4   0.1666667  0.1666667        0.3333333    0  
## Site.5   0.1666667  0.1666667        0.3333333    0  
## Site.6   0.0000000  0.5555556        0.5555556    +  
## Site.7   0.2222222  0.7777778        1.0000000    +  
## Site.8   0.2500000  0.2500000        0.5000000    0  
## Site.9   0.2000000  0.4000000        0.6000000    +  
## Site.10  0.0000000  0.3333333        0.3333333    +  
## Site.11  0.0000000  0.5000000        0.5000000    +  
## Site.12  0.4285714  0.0000000        0.4285714    -  
## Site.13  0.2222222  0.1111111        0.3333333    -  
## Site.14  0.5000000  0.5000000        1.0000000    0  
## Site.15  0.0000000  0.5000000        0.5000000    +  
## Site.16  0.2500000  0.2500000        0.5000000    0  
## Site.17  0.0000000  0.0000000        0.0000000    0  
## Site.18  0.2000000  0.0000000        0.2000000    -  
## Site.19  0.4000000  0.2000000        0.6000000    -  
## Site.20  0.1666667  0.1666667        0.3333333    0  
## Site.21  0.1111111  0.2222222        0.3333333    +  
## Site.22  0.0000000  0.2000000        0.2000000    +  
## Site.23  0.5000000  0.5000000        1.0000000    0  
## Site.24  0.2000000  0.4000000        0.6000000    +  
## Site.25  0.2857143  0.1428571        0.4285714    -  
## Site.26  0.0000000  0.0000000        0.0000000    0  
## Site.27  0.0000000  0.3333333        0.3333333    +  
## Site.28  0.0000000  0.6000000        0.6000000    +  
## 
## $BCD.summary
##  mean(B/den) mean(C/den)   mean(D)   B/(B+C)   C/(B+C) Change
##    0.2047902   0.2866497 0.4914399 0.4167147 0.5832853    +  
## 
## $t.test_B.C
##                 mean(C-B)     Stat   p.param p.perm   p<=0.05
## Paired t.test  0.08185941 1.437708 0.1620069  0.151          
## 
## $BC
## [1] NA
## 
## attr(,"class")
## [1] "TBI"
# Calculando Jaccard de forma independiente
betadiver(rbind(colSums(mc_marzo_para_tbi), colSums(mc_mayo_para_tbi)), method = 'j')
##       1
## 2 0.625
# Para contrastar con tabla de valores significativos.

Resumen de resultado sobre TBI y diversidad Beta estacional

Calculando distintos indices de disimilaridad (%diferencias, Jaccard, Sorensen), y aplicando pruebas de permutacion a la diversidad Beta temporal (T1 y T2), se obtiene un resultado consistente sobre un cambio positivo en la diversidad de especies (ganancia de especies), pero este cambio no es significativo en terminos estadisticos.

Asimismo, tras calcular el indice de Jaccard mediante betadiver, contrastamos con la tabla de valores significativos de Raymundo Real, evidenciamos que existen diferencias significativas en la diversidad Beta comparando las matrices combinadas (pooled) de las estaciones para un nivel de significancia de 0.05, pero no para un nivel de significancia de 0.01. Esto sugiere que las diferencias en cuanto a la composicion de hormigas nidificantes entre ambas estaciones son realmente bajas.

Analisis univariado de la precipitacion, meses marzo y mayo de 2022

precip_orig <- read_excel('JARDIN BOTANICO (JOSE RAMON).xls', range = "A9:I40")
precip <- precip_orig %>% 
  mutate_all(~ gsub('INAP', '0', .x)) %>%
  mutate_all(~ gsub('-', NA, .x)) %>% 
  mutate_all(~ as.numeric(.x))
sapply(precip, summary)
## $DIA
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     1.0     8.5    16.0    16.0    23.5    31.0 
## 
## $ENE
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.0000  0.0000  0.0000  0.2484  0.0000  3.3000 
## 
## $FEB
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   0.000   0.000   1.518   1.925  11.400       3 
## 
## $MAR
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   1.000   7.116  13.200  34.800 
## 
## $ABR
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   0.000   0.000   7.323   6.025  63.200       1 
## 
## $MAY
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   0.000   3.035   2.250  29.500 
## 
## $JUN
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   0.000   0.000   3.843   1.300  68.800       1 
## 
## $JUL
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   1.000   9.097  13.950  47.800 
## 
## $AGO
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   0.000   0.000   0.000   7.353  17.325  34.500       1
sapply(precip, sum, na.rm=T)
##   DIA   ENE   FEB   MAR   ABR   MAY   JUN   JUL   AGO 
## 496.0   7.7  42.5 220.6 219.7  94.1 115.3 282.0 220.6
precip_mar_colec <- precip %>% 
  select(DIA, MAR) %>% 
  filter(DIA %in% c(21:25, 30, 31)) %>% 
  pull(MAR)
sum(precip_mar_colec)
## [1] 27.1
precip_may_colec <- precip %>% 
  select(DIA, MAY) %>% 
  filter(DIA %in% c(10:13, 16:19, 31)) %>% 
  pull(MAY)
sum(precip_may_colec)
## [1] 46.7

Ecologa espacial

Generar hexagonos topologicamente correctos

Autocorrelación espacial mediante prueba Mantel (matrices de distancia)

library(sp)
# Matriz transformada Hellinger
mc_nidos_hell <- decostand(mc_nidos, 'hellinger')
# Matriz ambiental con centroides
st_geometry(hex_inters) <- 'geom'
ma_sf <- hex_inters %>% 
  inner_join(ma %>%
               rownames_to_column('ENLACE') %>%
               mutate(ENLACE = as.integer(ENLACE)))
ma_sp <- ma_sf %>% as_Spatial()
centroides <- ma_sf %>% st_centroid
ma_xy <- centroides %>% st_coordinates %>% as.data.frame
(vecindad <- centroides %>% dnearneigh(d1 = 0, d2 = 75))
## Neighbour list object:
## Number of regions: 28 
## Number of nonzero links: 128 
## Percentage nonzero weights: 16.32653 
## Average number of links: 4.571429
plot(ma_sp)
plot(vecindad, coords = ma_xy, add=T, col = 'red')

(pesos_b <- nb2listw(vecindad, style = 'B'))
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 28 
## Number of nonzero links: 128 
## Percentage nonzero weights: 16.32653 
## Average number of links: 4.571429 
## 
## Weights style: B 
## Weights constants summary:
##    n  nn  S0  S1   S2
## B 28 784 128 256 2544
# Matriz sin tendencia
mc_nidos_sin_tendencia <- resid(
  lm(as.matrix(mc_nidos_hell) ~ .,
     data = ma_xy))
(autocor_global_residuos_mc <- sapply(
  dimnames(mc_nidos_sin_tendencia)[[2]],
  function(x)
    moran.mc(
      x = mc_nidos_sin_tendencia[,x],
      listw = pesos_b,
      zero.policy = T,
      nsim = 9999),
    simplify = F))
## $`Pheidole subarmata`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.12312, observed rank = 2306, p-value = 0.7694
## alternative hypothesis: greater
## 
## 
## $`Dorymyrmex antillanus`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.084899, observed rank = 3532, p-value = 0.6468
## alternative hypothesis: greater
## 
## 
## $`Paratrechina longicornis`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.13807, observed rank = 1794, p-value = 0.8206
## alternative hypothesis: greater
## 
## 
## $`Solenopsis geminata`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = 0.14283, observed rank = 9320, p-value = 0.068
## alternative hypothesis: greater
## 
## 
## $`Pheidole jelskii`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.16881, observed rank = 1021, p-value = 0.8979
## alternative hypothesis: greater
## 
## 
## $`Pheidole moerens`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.093401, observed rank = 1420, p-value = 0.858
## alternative hypothesis: greater
## 
## 
## $`Mycetomoellerius jamaicensis`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.082011, observed rank = 3406, p-value = 0.6594
## alternative hypothesis: greater
## 
## 
## $`Tetramorium lanuginosum`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.093645, observed rank = 2933, p-value = 0.7067
## alternative hypothesis: greater
## 
## 
## $`Odontomachus ruginodis`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.29487, observed rank = 30, p-value = 0.997
## alternative hypothesis: greater
## 
## 
## $`Wasmannia auropunctata`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.093884, observed rank = 3161, p-value = 0.6839
## alternative hypothesis: greater
## 
## 
## $`Cyphomyrmex rimosus`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = 0.024807, observed rank = 7402, p-value = 0.2598
## alternative hypothesis: greater
## 
## 
## $`Tetramorium bicarinatum`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.097105, observed rank = 1337, p-value = 0.8663
## alternative hypothesis: greater
## 
## 
## $`Brachymyrmex heeri`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.096045, observed rank = 287, p-value = 0.9713
## alternative hypothesis: greater
## 
## 
## $`Pogonomyrmex schmitti`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = -0.19293, observed rank = 531, p-value = 0.9469
## alternative hypothesis: greater
## 
## 
## $`Monomorium floricola`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = 0.018724, observed rank = 8068, p-value = 0.1932
## alternative hypothesis: greater
## 
## 
## $`Acropyga dubitata`
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mc_nidos_sin_tendencia[, x] 
## weights: pesos_b  
## number of simulations + 1: 10000 
## 
## statistic = 0.018724, observed rank = 8091, p-value = 0.1909
## alternative hypothesis: greater
(autocor_global_residuos <- sapply(
  dimnames(mc_nidos_sin_tendencia)[[2]],
  function(x)
    moran.test(
      x = mc_nidos_sin_tendencia[,x],
      listw = pesos_b,
      zero.policy = T),
    simplify = F))
## $`Pheidole subarmata`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -0.77838, p-value = 0.7818
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.12312262       -0.03703704        0.01223152 
## 
## 
## $`Dorymyrmex antillanus`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -0.41818, p-value = 0.6621
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.08489939       -0.03703704        0.01309983 
## 
## 
## $`Paratrechina longicornis`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -0.92421, p-value = 0.8223
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.13807459       -0.03703704        0.01195157 
## 
## 
## $`Solenopsis geminata`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = 1.576, p-value = 0.05752
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.14282977       -0.03703704        0.01302614 
## 
## 
## $`Pheidole jelskii`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -1.1842, p-value = 0.8818
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.16880693       -0.03703704        0.01238092 
## 
## 
## $`Pheidole moerens`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -1.0731, p-value = 0.8584
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.093401039      -0.037037037       0.002759073 
## 
## 
## $`Mycetomoellerius jamaicensis`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -0.50945, p-value = 0.6948
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.082011111      -0.037037037       0.007793372 
## 
## 
## $`Tetramorium lanuginosum`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -0.62151, p-value = 0.7329
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.093645079      -0.037037037       0.008295785 
## 
## 
## $`Odontomachus ruginodis`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -2.3369, p-value = 0.9903
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.29487015       -0.03703704        0.01217280 
## 
## 
## $`Wasmannia auropunctata`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -0.50626, p-value = 0.6937
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.09388384       -0.03703704        0.01260874 
## 
## 
## $`Cyphomyrmex rimosus`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = 0.56388, p-value = 0.2864
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.02480714       -0.03703704        0.01202887 
## 
## 
## $`Tetramorium bicarinatum`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -1.1375, p-value = 0.8723
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##      -0.097105008      -0.037037037       0.002788338 
## 
## 
## $`Brachymyrmex heeri`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -1.904, p-value = 0.9715
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##     -0.0960451331     -0.0370370370      0.0009604943 
## 
## 
## $`Pogonomyrmex schmitti`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = -1.4674, p-value = 0.9289
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.19292906       -0.03703704        0.01128672 
## 
## 
## $`Monomorium floricola`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = 0.88658, p-value = 0.1877
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.018724330      -0.037037037       0.003955756 
## 
## 
## $`Acropyga dubitata`
## 
##  Moran I test under randomisation
## 
## data:  mc_nidos_sin_tendencia[, x]  
## weights: pesos_b    
## 
## Moran I statistic standard deviate = 0.88658, p-value = 0.1877
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.018724330      -0.037037037       0.003955756
pvalor_moran <- sapply(names(autocor_global_residuos),
                       function(x) autocor_global_residuos[[x]]$p.value)
tabla_patrones_espaciales <- data.frame(
  Especie = names(pvalor_moran),
  ValorP = pvalor_moran,
  Patron = ifelse(pvalor_moran<0.1, 'Agregado', 'Aleatorio'))
tabla_patrones_espaciales %>% kable
Especie ValorP Patron
Pheidole subarmata Pheidole subarmata 0.7818267 Aleatorio
Dorymyrmex antillanus Dorymyrmex antillanus 0.6620916 Aleatorio
Paratrechina longicornis Paratrechina longicornis 0.8223113 Aleatorio
Solenopsis geminata Solenopsis geminata 0.0575183 Agregado
Pheidole jelskii Pheidole jelskii 0.8818410 Aleatorio
Pheidole moerens Pheidole moerens 0.8583758 Aleatorio
Mycetomoellerius jamaicensis Mycetomoellerius jamaicensis 0.6947807 Aleatorio
Tetramorium lanuginosum Tetramorium lanuginosum 0.7328685 Aleatorio
Odontomachus ruginodis Odontomachus ruginodis 0.9902783 Aleatorio
Wasmannia auropunctata Wasmannia auropunctata 0.6936615 Aleatorio
Cyphomyrmex rimosus Cyphomyrmex rimosus 0.2864181 Aleatorio
Tetramorium bicarinatum Tetramorium bicarinatum 0.8723456 Aleatorio
Brachymyrmex heeri Brachymyrmex heeri 0.9715441 Aleatorio
Pogonomyrmex schmitti Pogonomyrmex schmitti 0.9288624 Aleatorio
Monomorium floricola Monomorium floricola 0.1876521 Aleatorio
Acropyga dubitata Acropyga dubitata 0.1876521 Aleatorio
# mc_nidos_sin_tendencia_d <- dist(mc_nidos_sin_tendencia)
# (mc_nidos_correlograma <- mantel.correlog(
#   mc_nidos_sin_tendencia_d,
#   XY = ma_xy,
#   nperm = 999))
# plot(mc_nidos_correlograma)